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Robustness of Majorana Modes and Minigaps in a Spin-Orbit-Coupled 
Semiconductor-Superconductor Heterostructure 

Li Mao and Chuanwei Zhangj 
Department of Physics and Astronomy, Washington State University, Pullman, WA, 99164 USA 

We study the robustness of Majorana zero energy modes and minigaps of quasiparticle excitations 
in a vortex by numerically solving Bogoliubov-deGennes equations in a heterostructure composed of 
an s-wave superconductor, a spin-orbit-coupled semiconductor thin film, and a magnetic insulator. 
This heterostructure was proposed recently as a platform for observing non-Abelian statistics and 
performing topological quantum computation. The dependence of the Majorana zero energy states 
and the minigaps on various physics parameters (Zeeman field, chemical potential, spin-orbit cou- 
pling strength) is characterized. We find the minigaps depend strongly on the spin-orbit coupling 
strength. In certain parameter region, the minigaps are linearly proportional to the s-wave super- 
conducting pairing gap As , which is very different from the A^ dependence in a regular s- or p-wave 
superconductor. We characterize the zero energy chiral edge state at the boundary and calculate 
the STM signal in the vortex core that shows a pronounced zero energy peak. We show that the 
Majorana zero energy states are robust in the presence of various types of impurities. We find the 
existence of impurity potential may increase the minigaps and thus benefit topological quantum 
computation. 

PACS numbers: 03.67.Lx, 71.10.Pm, 74.45.-l-c 



I. INTRODUCTION 

Topological quantum computation (TQC)^'^, where 
quantum information is processed using a decoherence- 
free subspace guaranteed by topological order, is a revo- 
lutionary new alternative to conventional quantum com- 
puting. In TQC, quantum information is encoded in 
certain nonlocal, topological, degrees of freedom of the 
underlying physical system (i.e., hardware) that do not 
couple to weak local noise. This special hardware, 
called 'non-Abelian topological matter', has been pro- 
posed to exist in certain classes of two-dimensional 
(2D) strongly-correlated systems, such as the v — 5/2 
fractional quantum Hall (FQH) systen>2ri^, chiral p- 
wave superconductor/superfluidiSr^, and some artificial 
states of cold atoms in optical lattices^i"— . In these sys- 
tems, the ground state wave function is a linear combina- 
tion of states from a degenerate subspace, and a pair- wise 
exchange of the particle coordinates unitarily rotates the 
ground state wave function in this subspace. Therefore, 
the exchange statistics of the particles is given by a multi- 
dimensional unitary matrix representation (as opposed 
to just a phase factor for bosons and fermions) of the 2D 
braid group, and the statistics is non-Abelian. 

Despite the tremendous technological potential, non- 
Abelian topological matter is rare in nature and gen- 
erally hard to observe in experiments'^. To circumvent 
this problem, there has been considerable interests re- 
cently for exploring the possibility of 'designing' non- 
Abelian topological order in the fertile laboratory of 
the cold atom systems^i^ and the regular solid state 
materials^^""— . Two important resources for the emer- 
gence of non-Abelian statistics in a correlated matter are 
(a) chirality of the constituent particles and (b) super- 
conducting order. In a composite system, these two basic 
ingredients of topological order may arise from two dif- 



ferent physical effects (for instance, spin-orbit coupling 
and s-wave superconductivity, respectively) to design a 
non-Abelian quantum state. This strategy allows us to 
use s-wave superconductors/superfluids, which are much 
more abundant in nature and much less sensitive to dis- 
order effects than their p-wave counterparts. 

One important recent progress along this direction was 
the proposed solid state heterostructure composed of an 
electron-doped semiconductor thin film, an s-wave su- 
perconductor, and a magnetic insulator as a non-Abelian 
platform for TQCS^. In this heterostructure, the super- 
conducting pairing order is induced on the semiconduc- 
tor thin film from the s-wave superconductor through 
the superconducting proximity effect, and the chirality 
of particles is supplied by the Rashba spin-orbit coupling 
in the semiconductor™. With a set of vortices in the het- 
erostructure, there exist one Majorana zero energy state 
of quasiparticle excitations in each vortex core. These 
zero energy states are the topological, degenerate, ground 
states following non-Abelian statistics and can be used to 
perform TQC. 

In addition to the existence of Majorana zero energy 
modes, another key ingredient for the physical implemen- 
tation of TQC is that the degenerate ground state sub- 
space must be separated from other non-topological ex- 
cited states by an energy gap, so that finite temperature 
cannot populate the excited states and ruin the topolog- 
ical properties of the system. Recently, it was shown^, 
by constructing a new type of Majorana operators that 
contain both zero energy and excited states, that the 
topological braiding statistics of Majorana fermions may 
still preserve even in the presence of non-topological exci- 
tations. However, the minigap is still important because 
the signal strengths of measurements will be reduced sig- 
nificantly when the temperature is comparable with the 
minigap where the number of non-topological excitations 



become significant. The magnitude of tfie energy gap di- 
rectly determines the operating temperature of the un- 
derlying physical system as a realistic TQC platform and 
is of critical importance. Furthermore, both Majorana 
zero energy states and the magnitude of the gap must be 
robust in the presence of various impurities. Although 
the existence of Majorana zero energy states in the vor- 
tex core in a heterostructure has been proved by analyz- 
ing the zero energy solution of the Bogoliubov-deGennes 
(BdG) equation, the magnitudes of the minimum energy 
gaps and their robustness against impurities have not 
been addressed. 

In this paper, we numerically solve the BdG equation 
for a vortex in the heterostructure and calculate the mini- 
mum energy gap (minigap) between the zero energy state 
and the first quasiparticle excited state in the vortex core. 
In the simulation, the proximity-induced s-wave super- 
conducting pairing gap Ag in the semiconductor is ob- 
tained from the self-consistent solution of the BdG equa- 
tions for a pure s-wave superconductor. The main results 
are summarized as follows: 

1) The full numerical simulation of the BdG equations 
confirms that the Majorana zero energy states exist only 
in the parameter region Vz > y/Jj^+A^, which agrees 
with previous results obtained from an approximate an- 
alytical approach. Here Vz is the perpendicular Zeeman 
field induced by the proximity contact with the magnetic 
insulator, fi is the chemical potential of the electron gas. 

2) The dependence of the minigap Eg on various 
parameters (14, /z, a) is characterized. Here a is the 
strength of the Rashba spin-orbit coupling in the semi- 
conductor. We find Eg depends strongly on a. In certain 
parameter region. Eg is ^ As , instead of ^ A^ for a reg- 
ular chiral p-wave superconductor/superfluid^. 

3) An analytical theory is developed to explain the 
properties of the zero energy chiral edge states at the 
boundary. 

4) The scanning tunneling microscopy (STM) signals 
around the vortex show pronounced peaks at the zero 
energy and the minigap, thus can be used to detect the 
zero energy modes and the minigaps in experiments. 

5) We show that Majorana zero energy modes are ro- 
bust in the presence of various impurity potentials. Sur- 
prisingly, we find the existence of impurity potential in 
the vortex core may increase the magnitudes of the mini- 
gaps, and thus may be useful for the physical implemen- 
tation of TQC. 

The paper is organized as follows: Section II lays 
out the BdG equation for a vortex in a semiconductor- 
superconductor heterostructure. Section HI discusses the 
parameter dependence of the zero energy states and the 
minigap Eg. In Sec. IV, we discuss the robustness of 
the zero energy states and the minigaps in the presence 
of realistic impurities. Section V consists of conclusions. 
The details about the numerical approach to the BdG 
equations are presented in Appendix A. The finite size 
effect in the BdG equations is discussed in Appendix B. 
In Appendix C, we discuss the zero energy chiral edge 
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Figure 1. (Color online) (a) An illustration of the structure 
of the spin-orbit coupled semiconductor-superconductor het- 
erostructure. (b) The single particle energy spectrum in a 
spin-orbit coupled semiconductor. 



modes at the boundary. 



II. BDG EQUATIONS FOR A VORTEX 

The physical system we consider is a heterostructure 
composed of an s-wave superconductor, an electron- 
doped semiconductor thin film, and a magnetic insulator 
(Fig. [T^). The dynamics of electrons in the semiconduc- 
tor are described by a single particle effective Hamilto- 
nian 



Hn = 



P 

2m' 



- a{px(Ty-Py(Jx)+Vz(Tz- ti, (1) 



where m* is the conduction-band effective mass of elec- 
trons, /i is the chemical potential, a is the strength of the 
Rashba spin-orbit coupling, Vz is a perpendicular Zeeman 
field induced by the proximity contact with the magnetic 
insulator, ai are the Pauli matrices for the electron spins. 
The Hamiltonian yields two spin-orbit bands (Fig. [TJd) 
with energy dispersions 
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in a uniform system. Henceforth we set h = \. A fi- 
nite energy gap TVz is opened at fc = for a nonzero Vz- 
When the chemical potential lays in the gap, electrons 
only occupy the lower spin-orbit band at a low tempera- 
ture. 

The mean field Hamiltonian of an s-wave superconduc- 
tor can be written as 



U 



dx 



( 4 ^^^ "^ ("^ 



H, A(r) 
A*(r) ~H, 



at (r) 
4(r) 



(3) 

in the Nambu space, where ils = — V'^/2r7i — ii^i? + C^(r) is 
the single particle Hamiltonian, Ep is Fermi energy, C/(r) 
is an external potential, a^ (r) are annihilation operators 
of electrons for position eigenfunctions rather than for 
momentum eigenfunctions, A(r) is the s-wave pairing or- 
der parameter. The Hamiltonian ([3]) can be diagonalized 



by the Bogoliubov transformation 



\{v))-^{Z{r) K(v) }[^lj^ (4) 



where the wavefunctions u„(r), ti„(r) satisfy the BdG 
equation 



Hs A(r) \ (un{v)\ ^ (un{v) 
A*(r) -hJ {v„{r)J '\v„{r) 



(5) 



and the normahzation condition 

dr[<j(r)u„(r) + w*,(r)u„(r)] = ^„j„. (6) 

The order parameter A(r) can be determined 
self-consistently^ through the relation A(r) — 
5^„u„(r)w*(r). Here g is the effective electron- 
electron interaction strength in the superconductor. 

The proximity effect between the s-wave superconduc- 
tor and the semiconductor induces an effective supercon- 
ducting pairing for electrons in the semiconductor de- 
scribed by the Hamiltonian 



Hp^ f dr { A,(r)4 (r) c| (r) + H.c.} 



(7) 



where cj. (r) are the creation operators for electrons, and 
As(r) is the proximity- induced effective s-wave pairing 
gap in the semiconductor thin film. Because of the 
Rashba spin-orbit coupling in the semiconductor, the 
single particle Hamiltonian Hg in the BdG equation ([5]) 
should now be replaced with Hq in ([T]). The BdG equa- 
tion written in the Nambu spinor basis becomes 



Ho 

A:(r) 



A,(r) 



$„(r) - £;„$„(r) 



(8) 



where $„(r) = ['u„|(r), u„;(r),'y,4(r), -u„|(r)] is the 
quasiparticle wavefunction. The Bogoliubov quasiparti- 
cle operator is 



7n = h^^Yl [Una{r)ci (r) + Vna{r)Ca (r)] 



(9) 



In the presence of a vortex in the semiconductor- 
superconductor heterostructure, the order parameter 
takes the form As{r,d) — As(r)e'^, and the solutions 
(£'„,$„(r)) of the BdG equation ^ correspond to the 
quasiparticle excitation energies and states in the vor- 
tex core. For simplicity of the calculation, we consider 
a two-dimensional cylinder geometry with a hard wall at 
the radius r ~ R and a single vortex at r = 0. This 
system preserves the rotation symmetry and the BdG 
equation can be decoupled into different angular momen- 
tum channels indexed by I with the corresponding spinor 
wavefunction 



<(r) = e^'« [uU{r),ul,ir)e^o,vl,{r)e-^<^,-vl^{r)]- 




(10) 



Figure 2. Plot of the s-wave pairing gap As (r) from the self- 
consistent solution of the BdG equation ((5]) for a pure s-wave 
superconductor. 



Note that the BdG equation has the particle-hole sym- 
metry, therefore if $Ji(r) is a solution with an energy 
E, then there is another solution with the energy —E in 
the angular momentum — / channel. Henceforth we only 
consider E > solutions. 



III. MAJORANA MODES AND MINIGAPS 

We numerically solve the BdG equation ([8]) with a vor- 
tex and calculate the quasiparticle excitation energy _E„ 
and wavefunction $,i(r) for various parameters {Vz, fJ-, a). 
In the numerical treatment, the radial wave functions 
Mjjg.(r), wJi(j(r) in Eq. (|10l) are expanded on an orthogonal 
basis (pji{r) = y/2Ji{l3jir/R)/[RJi+i{(3ji)], where J; (a;) is 
the l-th order Bessel function, (3ji is the j-th zero of Ji{x). 
This basis satisfies the boundary condition 4)ji{R) — 
automatically. The BdG Hamiltonian ([S]) can be writ- 
ten as a matrix form on this basis and then diagonalized 
to obtain En and $n(r) (more details about the numeri- 
cal method can be found in Appendix A). Henceforth we 
choose fc^^ — kp^ (m/m*) ' — 5kp^ as the length unit 
and 77 — h'^k'^/2m* = Ep as the energy unit, where kp is 
the Fermi wavevector in the s-wave superconductor, and 
we adopt an effective mass m* = 0.04771 for electrons in 
the conduction band of a semiconductor. 

The s-wave pairing gap A5(r) in the BdG equation 
(|S]) for the semiconductor-superconductor heterostruc- 
ture is obtained by solving the BdG equation ^ self- 
consistently for a pure s-wave superconductor—, and the 
resulting A^ (r) is plotted in Fig. [2j In the self-consistent 
procedure, we first guess a form of As(r) and inset it 
into the BdG Eq. (O, from which we obtain the wave- 
function Un{r) and Vn{r). As(r) is then self-consistently 
determined from u„(r) and Vn{r) through the relation 
As(r) — ffX^Ti '"n(^)^n(^)- ^^^ '^^^ As(r) is inset into 
the BdG Eq. ([5]) to start another cycle of the calcu- 
lation. The procedure continues until As(r) converges. 
As(r) is zero at the vortex core, approaches a uniform 
value Ao « 0.11 in the bulk, and drops to zero at the 



boundary, as expected. 

The finite size effect can become important in tlie cal- 
culation of the Majorana zero energy modes and the mini- 
gaps in certain parameter region in the numerical simu- 
lation of the BdG equation. More details about the finite 
size effect can be found in Appendix B. Here we choose a 
large radius R = 250 of the cylinder to suppress the finite 
size effect. In practice, solving the BdG equation self- 
consistently to obtain A^ (r) for a pure s-wave supercon- 
ductor with a large radius R = 250 is very time costing. 
Since Ag (r) approaches the uniform bulk value Ao within 
a finite distance from the vortex core, we choose a pair- 
ing gap As(r) based on the numerical result As(r) from 
the R = 5k~^ calculation (Fig. O, i.e., As{r) = As{r) 
for < r < 2.8, As(r) = Aq for 2.8 <r < R~ 2.2, and 
A,(r) ^As{r-R + 5) iov R- 2.2 < r < R. 

In Fig. [3l we plot the quasiparticle energy Eni at dif- 
ferent angular momentum I channels. We see only at the 
I = channel, there is a unique Majorana zero energy 
solution. At non-zero I channels, there are two discrete 
energy levels: one is the edge state, the other is the vor- 
tex core state. This can be clearly seen from the corre- 
sponding eigenwavefunctions for these two energy levels 
at the I = —1 channel, which are plotted in Fig. |4l The 
continuous spectrum in Fig. |3] corresponds to the bulk 
excitations. Inside the vortex core, the first excited state 
above the zero energy mode is at the I = —1 channel and 
the corresponding energy difference is called the minigap. 
The minigap is the energy gap that protect the Majorana 
zero energy state from finite temperature and disorder ef- 
fects, therefore its magnitude is crucially important for 
the observation of non- Abclian statistics of quasiparticles 
in this system. 

The non-Abelian topological property of the Majorana 
zero energy states at the I = channel originates from 
their special forms of the wavefunctions Ua-{r), Va{r), 
which lead to a self-Hcrmitian Bogoliubov quasiparticle 
operator, 7q = 70. Specifically, the quasiparticle wave- 
functions around the vortex core satisfy Ucr{r) = —Va{r), 
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Figure 3. (Color online) Plot of the quasiparticle energies E„t 
in a vortex for the angular momentum I from —4 to 4. a — 1, 
11 = 0, and K = 0.3. 
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Figure 4. Plot of the wavefunctions v-f (r) of the vortex core 
state and the edge state in the I — —1 angular momentum 
channel, fj. = 0, Vz = 0.3, a = 1. 
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Figure 5. (Color online) Plot of the wavefunctions Uc (r), 
Va (r) of the Majorana zero energy state, fi = 0, Vz = 1, 
a = 1. The wavefunctions have two parts. Around r — 
is the zero energy state in the vortex core. Around the edge 
r = _R is the zero energy chiral edge state. In the vortex core, 
u^ (r) = - v^ (r). At the edge, w^ (r) = v^ (r). 



as clearly seen from Fig. \5\ Chosen an artificial over- 
all phase e*'^/^ for the wavefunction ([TU]) . it is easy to 
show the Bogoliubov quasiparticle operator 70 defined 
in Eq. (jH]) satisfies 7q = 79. Note that 70 cannot be 
the electron creation or annihilation operator because it 
does not obey the anticommutation relation for electrons. 
Instead , the quasiparticle defined by 70 is called a Ma- 
jorana fermion. The exchange statistics of two Majorana 
quasiparticle excitations in two vortex cores are not the 
simple Fermi statistics: they can be either Abelian or 
non-Abelian in the degenerate subspace spanned by the 
Majorana operators in a set of vortices^. 
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Figure 6. (Color online) Plot of three quasiparticle excitation 
energies En with respect to the Zeeman field Vz in a vortex 
core. Solid line: ground state energy at Z = 0. Dashed line: 
ground state energy at Z = — f . Dotted line: the first excited 
state energy at Z = 0. a = 1, /x = 0.2. 



Interestingly, the wavefunctions still oscillate spatially 
even deep inside the bulk region. Such an oscillation 
makes it difficult to realize a single qubit gate for uni- 
versal TQC using the tunneling between two vortices^^. 
We also observe the zero energy edge states around the 
boundary, which satisfy Ucir) = Vc(r\ in contrast to 
u^{r) — —v^{r) in the vortex core. The wavefunctions 
Ua-{r) and Va{r) vanish at the boundary, as expected. An 
analytic explanation for the observed u^(r) = Va{r) re- 
lation of the chiral edge modes is provided in Appendix 
C. We emphasize that zero energy state in the vortex 
core and the edge must appear simultaneously because 
Majorana zero energy modes only come in pairs (either 
between two vortices or between a vortex and the edge). 

In Figs. (|6l7l8p . we plot three different quasiparticle 
energies with respect to various physical parameters: (i) 
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Figure 8. (Color online) Plot of the bulk excitation gap 
(dashed line) and the minigap (solid line) with respect to the 
spin-orbit coupling strength a. /^ = 0, 14 = 0.3. 



the ground state energy in the vortex core at the / = 
channel. In certain parameter region, this state is the 
zero energy Majorana mode; (ii) the ground state energy 
in the vortex core at the I = —1 channel. In the pa- 
rameter region with the zero energy modes, this energy 
corresponds to the minigap Eg] (iii) the first excited en- 
ergy at the I = channel. As we can see from Fig. [31 it 
corresponds to the minimum bulk excitation energy. Fig. 
ini shows the dependence of these quasiparticle energies on 
the Zeeman field V^. We see the Majo rana zero energy 
state exists only in the region Vz > -y^Ap + jj? k, 0.23. 
In this region, the minigap Eg first increases rapidly to a 
level E^'^^, and then decreases slowly with increasing Vz- 
With a very large Vz , all electrons occupy the spin down 
states, and there is no superconducting pairing. There- 
fore the decrease of the minigap with increasing Vz is ex- 
pected. Note that Eg is larger than the typical minigap 



Figure 7. (Color online) Plot of three quasiparticle excitation 
energies En with respect to the chemical potential /^ in a 
vortex core. The notation for different lines are the same as 
that in Fig. a = 1, K = 0.3. 



kQ ^ 0.01 for a regular (s-wave or chiral p-wave) 
superconductor/superfluid^. To characterize the depen- 
dence of the minigap Eg on the uniform bulk pairing gap 
Ao, we numerically calculate Eg for different Aq and plot 
it in Fig. El We see in the region Aq = 0.07 ~ 0.22, Eg is 
roughly proportional to Aq. When Aq > 0.22, the mini- 
gap decreases because Aq is now very close to Vz — 0.3. 
In general, Eg may be a linear combination of Ag and 

In Fig. [71 we see the zero energy modes disappear in 
the region ^ > \/Vz — Aq. As expected, the minigap Eg 
has a maximum at fj, — 0. Eg only changes slightly when 
the chemical potential varies. In Fig. El we plot the 
quasiparticle energy En with respect to the spin-orbit 
coupling strength a. We see the minigap has a strong 
dependence on a. Eg initially increases quickly with a 
growing a, and reaches the maximum, then decreases 
very slowly for large a. This is expectable because the 
Rashba spin-orbit coupling provides the necessary chi- 
rality for the zero energy states. When a — 0, the cou- 
pling between spin up and down states vanishes, and a 
pure s-wave superconductor does not have the zero en- 



ergy modes. 

The zero energy modes in the vortex core can be 
probed by bringing a STM tip close to the vortex core in 
experiments. The resuhing tunnehng conductance can 
be written as^ 



^ ^-Y,[uUr)nE,-eV)+vl{r)f'{E,+eV)], (11) 



where eV is the voltage bias, i represents different energy 
level, / is the tunneling current, T is the operation tem- 
perature of the STM, / = l/(exp((^, - eV) /ksT) + 1) 
is the Fermi-Dirac distribution, and the derivative of / is 
with respect to E. In Fig. [TUl we plot the STM tunneling 
conductances at different radius of the vortex core, which 
show clear zero energy peaks coming from the zero energy 
states. The other peaks correspond to other vortex core 
states at the angular momentum I ^ channels. The first 
peaks around the zero energy are at the I — ~1 channels. 
The distance between the peak centers at I = 0, —1 chan- 
nels is a measurement of the minigap. There is also an 
asymmetry of the peak strengths at zLeV. Because of the 
symmetry E — >■ —E, uiaij) o V-iair) in the BdG equa- 
tion, we only need the Ei > terms in Eq. ([Tl]) . There- 
fore the magnitudes of the peaks at positive or negative 
eV are proportional to X^o- ""■fa- ('') ^^'^ J^a'^ia-i''") respec- 
tively. Generally, J^a'^fai''') ^ Y.a'^lA'r), leading to the 
asymmetric peaks around eV = 0. In particular, we find 
the peak for r = at the negative eV disappears, which 
means ^^ i'L(O) must be zero for the minigap state (i.e., 
Z = — 1, E^i > state). As we can see from Eq. (|Aip . 
X^ii|^(r) in the minigap state involves the basis states 



4>i=~2,j {r) and (pi^^i.j (r) which are zero at r 
cause the Bessel function J/^o (^ = 0) = 0. 
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Figure 10. (Color online) Plot of the STM tunneling conduc- 
tance at different r around the vortex, a = 1, /i = 0, 14 = 0.3, 
ksT = Ao/50. 
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Figure 11. (Color online) Plot of three quasiparticle excitation 
energies En with respect to the Zeeman field 14 in the vortex 
core in the presence of a Gaussian impurity. The notation for 
different lines are the same as that in Fig. [6] /i = 0.2, a = 1. 
Uo = 10, s = 0.5. 
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Figure 9. (Color online) Plot of the minigap (solid line) and 
the bulk excitation gap (dashed line) with respect to the uni- 
form bulk superconducting pairing gap Aq. /i = 0, 14 = 0.3, 
a = \. 



IV. EFFECTS OF IMPURITIES 

In a realistic TQC platform, both Majorana zero en- 
ergy states and minigaps need be robust in the presence 
of impurities. A general argument on the robustness of 
zero energy states is based on the particle-hole symmetry 
in the BdG equation, which ensures its energy spectrum 
to be symmetric upon E -^ —E, that is, the E and —E 
states must come in pairs. However, the zero energy 
does not come in pairs in a BdG equation (see Fig. [3]). 
Therefore a local small perturbation cannot destroy this 
symmetry and shifts the zero energy to a finite energy E 
that must emerge simultaneously with another state with 
an energy —E. Currently, the robustness of the magni- 
tude of the minigap in the presence of impurities is still 
not clear. In the following, we consider three types of 
rotationally invariant impurities and study their effects 
on the zero energy states and the minigaps in the vortex 
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Figure 12. (Color online) Plot of the bulk excitation gap 
(dashed line) and the minigap (solid line) with respect to the 
spin-orbit coupling strength a (a) and the impurity strength 
Uo (b) in the presence of a Gaussian impurity potential. 14 = 
0.3, ^ = 0, s = 0.5. (a) f/o = 10. (b) a = 1. 



core. Interestingly, we find that impurities may actually 
increase the minigaps in the vortex. 

First, we consider a spin- independent Gaussian impu- 
rity 



U{r) = C/o exp (-rV2s^) 



(12) 



in the BdG equation ([8]) for a vortex, where C/o is the 
impurity strength, s is the half- width of the impurity po- 
tential that is comparable to the size of the vortex core. 
In Fig. [Til 'W6 plot the quasiparticle energies with re- 
spect to Vz in the presence of a Gaussian impurity (|12p . 
In a regular s-wave superconductor, it is expected that 
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Figure 13. (Color online) Plot of three quasiparticle excitation 
energies E„ with respect to the Zeeman field 14 in the vortex 
core in the presence of a magnetic Gaussian impurity. The 
notation for different lines are the same as that in Fig. (6] 
/i = 0.2, a = 1, s = 0.5. (a) Uo = 10; (b) f/o = -10. 
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Figure 14. (Color online) Plot of three quasiparticle excitation 
energies En with respect to the Zeeman field 14 in the vortex 
core in the presence of a Gaussian random impurity (a) and 
a magnetic Gaussian random impurity (b). The notation for 
different lines are the same as that in Fig. [5] /^ = 0.2, a = 1. 
C/o = 1, s = 0.5. 



such an impurity potential can couple different quasi- 
particle excitation states, and thus modify the energy 
spectrum. This can be seen in Fig. [TT] in the param- 
eter region V^ < y/n'^ + Aq without the zero energy 
states. We find that the zero energy states in the region 
Vz > -y^/i^ -I- Aq are very robust even for a large impurity 
potential Uq = 10. More interestingly, we find that, by 
comparing with Fig. [51 the presence of the Gaussian im- 
purity potential can increase the magnitude of the mini- 
gaps significantly. In Fig. 1121 we plot the minigap and 
the bulk excitation gap in the presence of a Gaussian im- 
purity potential. We see the minigap approaches the bulk 
excitation gap for a large Uo- This enhancement of the 
minigap can be understood by considering the fact that 
the impurity potential can repulse (or attract) electrons 
and enlarge the energy level splitting between different 
discrete states in the vortex core. In practice, we may 
add a Gaussian type of potential at the vortex center 
to increase the minigaps and ensure the corresponding 
topological protection from finite temperature effects. 

Secondly, we consider a magnetic impurity potential 
U{r)az , whose effects on the zero energy modes and mini- 
gaps are shown in Fig. 1131 Such an impurity potential 
acts as a spatially localized Zeeman field. No matter the 
local Zeeman field U{r)az is along the same (Fig. [T3k ) 
or opposite (Fig. 113b ) direction as the original Zeeman 
field Vz, the zero energy states do not change and the 
minigaps are enhanced significantly. Note that in both 
cases the impurity does have a significant i mpact on the 
quasiparticle energies in the region Vz < y/n^ + Ag. 

Finally, we consider random impurity potentials 



^ ran v j 



JV 

E 



t/osinWmrexp(-r^/2s^) (13) 



and Uran{''')(Jz, where uim are random frequencies (see 
the inset in Fig. [Tib for the spatial profile of the random 
potential) . This potential still has significant amplitudes 
outside the vortex core region. In Fig. [TH we find that 
the critical T4 for the existence o f the zero energy modes 
is smaller than the expected y/JJ^+A^. It may orig- 
inate from the long range property of the random po- 
tential, whose average U around the vortex core shifts 
the chemical potential fj,. Therefore the critical point is 
shifted to ^/(/z — C/)^ -I- Aq. For a magnetic random im- 
purity Uran{'i')(Jz, the Corresponding critical point shifts 
to \//?+A^ — tj . However, the minigaps are small in 
these parameter regions, therefore they cannot be used as 
a TQC platform eve n the zer o energy states exist. Only 
in the region Vz > \/ ^^ -I- Aq, the zero energy states and 
the minigaps are robust against impurities. 

So far the impurity potentials have been chosen to be 
rotationally invariant to simplify the numerical calcula- 
tion. Non-rotationally invariant potentials couple differ- 
ent angular momentum channels and thus are difficult to 
simulate numerically. In a realistic situation, vortices in 
superconductors are often pinned by impurities so that 
the centers of the impurities and vortices are the same to 
minimize the free energy of the system. Such overlapping 
of the centers provides one justification of our choice of 
rotationally invariant impurity potentials. More gener- 
ally, an impurity potential may be expanded as 



V{r,Q)^ J2 C/„(r)exp(zn0), 



(14) 



where UQ{r) is the rotationally invariant part, while 
the Un^olr) corresponds to the non-rotationally in- 
variant part. For instance, a Gaussian impurity po- 
tential located at (a;oiO) instead of (0,0) can be 

written as U{r) = 

2 I 2 



exp 






exp 



2s^ 



exp I 



■ COS! 



il^ ) A^n=-oo-'ni^)(i^Piin9), where /„ is the 
modified Bessel functions. Because the eigenstates of 
the BdG equation ([8]) have fixed angular momentum /, 
Uo{r) only couples states with the same I, while Un^oif) 
couples states with different I. We can treat the im- 
purity potential as a perturbation to the original BdG 
Hamiltonian. The first order energy correction to the 
zero energy state is AE^^^ = J2a i'^Ocr\U{r,6)\uoa) - 
{voa-\U{r,9)\vorT) — because mqo- — -^voa around 



vortex. 

|2 



The second order correction AE, 



(2) _ 



the 

— J2i \Uoi\^ /Ei = because of the coexistence of Ei and 
—Ei in the energy spectrum. Therefore the zero energy 
states are preserved by the particle hole symmetry in 
the BdG equation, as discussed at the beginning of this 
section. As to the minigap state, we see the first or- 
der correction A^J^^ = ($5^="^(r)| C/o(r) |$'="^(r)) only 
depends on the rotationally invariant part. The non- 
rotationally invariant part Un^oif) only yields a second 
or higher order correction to the minigap, which is gen- 
erally small and may be neglected. 



CONCLUSION 



In conclusion, we show that Majorana zero energy 
states and minigaps of quasiparticle excitations in a vor- 
tex are robust in a heterostructure composed of an s-wave 
superconductor, a spin-orbit-coupled semiconductor thin 
film, and a magnetic insulator. By numerically solving 
BdG equations with a vortex in a cylinder geometry, the 
dependence of Majorana zero energy states and minigaps 
on various physics parameters of the system (Zeeman 
field, chemical potential, spin-orbit coupling strength) is 
characterized. In certain parameter region, the minigap 
is proportional to the pairing gap As, instead of ~ A^ 
for a regular chiral p-wave superconductor/superfiuid. 
The existence of the zero energy chiral edge state at the 
boundary is demonstrated both analytically and numeri- 
cally. The STM tunneling conductance in the vortex core 
shows pronounced peaks at the zero energy as well as the 
minigap energy, which can be used to measure the zero 
energy state and the minigap. We show that the Majo- 
rana zero energy states are robust in the presence of var- 
ious types of impurities. Surprisingly, we find that impu- 
rity potentials may greatly enhance the magnitudes of the 
minigaps. Therefore they can be induced intentionally in 
the heterostructure to enhance the topological protection 
of the non-Abelian platform. We believe our character- 
ization of the behaviors of the zero energy modes and 
the minigaps in different parameter regions will help the 
design of a realistic semiconductor-superconductor het- 
erostructure system for the experimental observation of 
non-Abelian statistics and the physical implementation 
of TQC in the future. 
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Appendix A: Numerical method for solving the BdG 
equation 



We expand the radial wavefunctions ul^^{r),vl^^{r) 
in Eq. pUj) using the basis states (t>ij{r) — 

V2Ji{l3ijr/R)/[RJi+i{(3ij)], where /3y is the j-th zero of 
the Bessel function Ji{x). The BdG eigenvalue equation 
(|S]) reduces to a block diagonal matrix 
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Figure 15. (Color online) Plot of the first two quasiparticle 
energies at the I = channel with respect to the Zeman field 
14. M = 0, a = 1. (a) R = 25; (b) R = 100. 



where 



,Pl 



J± 



{T^h = (-^±y^-^^)5^,, 



{Si)ij = a r(j>u{r){dr + 



l + l. 



(A2) 
k+ij{r)dr, (A3) 



(A 



i)i] 



rAs (r) (t)u{r)(t)i^ij{r)dr, 



7/' - \ii^ • • • 7/' • • -F 

"rwr ~ i^nal^ t "'71a j J J ) 

7/ = [7/ • • • 7;' • • -F 



(A4) 

(A5) 
(A6) 



Appendix B: Finite size effect 

In an infinite large system, the zero energy modes al- 
ways exist in the parameter region Vz > y/^ + JJ^^. In 
Ref.— , the existence of the zero energy modes was proven 
by solving the BdG equation separately in and out the 
vortex core edge and matching the boundary condition at 
the vortex edge. A step function for the pairing gap has 
been used in this case to simplify the calculation. The 
zero energy modes exist when the number of unknown 
coefficients of the wavefunctions is equal to the number 
of independent constraints (vortex edge boundary con- 
ditions and normalization of the wavefunction) , which 
occurs in the parameter region Vz > \/A'^ + jj?. How- 
ever, there will be additional constraints in the system 
boundary, therefore zero energy solutions do not exist 
generally in a finite size system, that is, the quantum 
well state may destroy the zero energy modes when the 
size of the material or the confining potential is small 
enough. 

To describe the finite size effect, we consider the pa- 
rameter region Vz — >■ oo, where the spin-orbit coupling 
and proximity induced superconductivity are not impor- 
tant. In this region, the BdG matrix (|A1[) is diagonal. 



and the zero energy solution exists when 



(Tf). 



i?2 



14=0. 



(Bl) 



Here we take /i = 0. We see there are zero energy modes 
only at some special Vz = iSf^/R^. When Vz exceeds 
PIj^/R^, the lowest two eigenvalues are Vz — Pl^/R^ and 
Pin+i/R^ — Vz, and the lowest energy branch reaches its 
local maximum at 



'in^ 



i?2 



Vz^Vz 



i?2 



where the integer n is determined by 



0L 

i?2 



<Vz< 



Pln+l 
i?2 ■ 



(B2) 



(B3) 



Pin is the n-ih zero of the Bessel function Ji (a:;) and 
satisfies Pin+i ~ Pin + TT, therefore the maximum energy 

is 



^m.n.x — 



1^171+1 - f^ui _ ti'n/T^ 



2i?2 



R 



(B4) 



As an example, we consider the parameter a = \, 
/i = 0, Z = and use the pairing gap from the s-wave 
superconductor. In Fig. (fTSj) . we plot the lowest two 
quasiparticle energies with respect to Vz for two different 
sizes of the system R — 25, 100. We see when Vz is big 
enough, the fitting function [Emax = ''^^/V'z/R) is a good 
approximation to the oscillation amplitude of the ground 
state. While at a small V^, the amplitude is much smaller 
than TT\/V/ R. It is clear when R is big, the oscillation 
amplitude is strongly suppressed. In practice, we choose 
the parameters R = 250, Vz < 2, and the finite size ef- 
fects can barely be seen and thus be neglected, as clearly 
demonstrated in Fig. ^. 



Appendix C: Chiral edge states at the / = channel 

At the I = channel, the BdG equation for the zero 
energy state can be written as 






where 



° " ' "^ F{r) + ^-Vz- /i 



-adr 



(CI) 



(C2) 



r (r) = —77 (92 -I- ^dr). The BdG equation can be fur- 
ther reduced to a 2 x 2 matrix differential equation 



F (r) + Vz 
— A As (r) — ad, 






(C3) 
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where the parameter A = ±1 and iayTy^oii^) ~ i^^o{f) 
for nondegenerate zero energy states. These conditions 
yield ^o(r) — [uf,ui] and M(j(r) — Xvcr{r). It has 
been shown in Re f.'^" that only in the parameter region 
X = —1, Vz > y^Aq + /x^, there is an unique zero en- 
ergy solution. Here we show that a zero energy chi- 
ral edge state exists in the parameter region A = 1, 

As we can see from Fig. [21 the pairing gap As{r) 
decreases from its bulk value Ag to zero at the bound- 
ary. For simplicity, we approximate the radial depen- 
dence As(r) around the boundary with a step function: 
A (r) = for i? > r > i? - (5, and A (r) = Aq for 
r < R — S, where 6 is the coherence length of As(r) 
around the edge. Because R is very large, the BdG equa- 
tion (IC3I) reduces to 



*o(r) = (C4) 



-77^2 + v; - /i AA, (r) + adr 
-AAs (r) - adr -ridl -Vz- ^J. 



around the boundary. In the region R > r > R 
6, the solution of Eq. (IC4p is given by ^o(r) 
[u^, Ui] exp (z (r — R)) with the constraint 



-?7z2 + Vz- ^J. 



za 



-■qz^ -Vz- fJ. 



u^ 
Ui 



0. (C5) 



The characteristic equation for z is (772^ + MJ ~ 
V^ — z^o? = 0. There are four solutions of Eq. 
(|C4p which are well behaved at the edge: <l>i(j) = 



exp{zi{r — R)) {i = 1,2,3,4), where zi, Z2, 

Z3, Z4, are four solutions of Eq. (IC5[) . The full wave- 
function in the region R > r > R — 6 can be written as 
*o(?') = ci0i(r) -I- C2(t>2{r) + C303(r) -f CA(t>A{r). 

Far from the boundary, where A (r) = Ap , wc can 
expand the solution as a series in j^z- 



*oW 



exp {iz [R- r)) 



E 



(i?-rr 



(C6) 



where a„ are the corresponding spinors. The zeroth order 
coefficient oq satisfies the following equation: 



r]z' + Vz 

— AAq -I- iza rjz 



/i AAq — iza 

2 V 



M 



flQ = 0. 



(C7) 



The higher order coefficients a„ can be calculated from 
flo using a set of recursion relations. The characteristic 
equation has 4 complex root for z, same as that in Ref.— . 
Because Im[z„] < is required for a physical solution, 
there are three indepe ndent ro ots only in the parameter 
region A = 1, 14 > •\/Aq -|- /i^, which yields three inde- 
pendent coefficients. Together with the four independent 
coefficients ci, C2, C3, C4, in the region R> r > R — 6, and 
seven constraints (match of ^oii^) and ^o(^) at r = R—6, 
the boundary condition ^o(^) = 0, and the normaliza- 
tion of the wavefunction) , we can obtain a unique zero 
energy edge state at the boundary. However, the chiral 
edge state wavefunction satisfies Ua{r) ~ Va{r) instead of 
Uaii") = —Va{r) for the vortex core state, as clearly seen 
from Fig. [5l 
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